\(~\)
This notebook helps visualize and explore surface and groundwater observations stored in the National Water Information System (NWIS). It was written to perform automated pulls and interactive mapping and plotting of data as it was being collected and analyzed for a specific project, but it can be modified for any region where there is data in NWIS.
\(~\)
# watershed screening sites
screensites <- whatNWISsites(bBox = c(-70.448, 41.626, -70.392, 41.677))
nums <- "504|505|506|507"
screensites <- screensites %>%
filter(grepl(nums, station_nm))
# local demonstration sites
sites <- whatNWISsites(bBox = c(-70.405, 41.669, -70.395, 41.677))
# site info
baseurl1 <- "https://nwis.waterdata.usgs.gov/usa/nwis/qwdata/?site_no="
baseurl2 <- "&agency_cd=USGS"
siteinfo <- readNWISsite(sites$site_no) %>%
dplyr::select(site_no, station_nm, site_tp_cd, dec_lat_va, dec_long_va, dec_coord_datum_cd, alt_va, alt_datum_cd, well_depth_va)%>%
mutate(id = substr(station_nm, 9, 11)) %>%
mutate(weblink = paste(baseurl1,site_no,baseurl2,sep = " "))
# reclassify site types for mapping
siteinfo$site_tp_cd[grep("M01", siteinfo$station_nm)] <- "multilevel sampler (MLS)"
clust <- "505-0|517-0|524-0|525-0|526-0"
siteinfo$site_tp_cd[grepl(clust, siteinfo$station_nm)] <- "well cluster"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="GW")] <- "water table well"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="LK")] <- "surface water (pond)"
# add a treatment variable
siteinfo$treatment <- NA
controls <- c(505,508,514,515,520,521)
siteinfo$treatment[which(siteinfo$id %in% controls)] <- "control"
impacts <- c(517,522,523,524,525,526)
siteinfo$treatment[which(siteinfo$id %in% impacts)] <- "impact"
siteinfo$treatment[siteinfo$id == 509] <- "other"
# site numbers within 10ft of water table (as of Sept, 2021)
sn <- c("414010070235201", "414015070235001", "414015070235002", "414015070235301", "414015070235303", "414016070235001", "414016070235002", "414016070235101", "414016070235102", "414016070235103", "414017070235004", "414017070235005", "414017070235006", "414017070235008", "414017070235010", "414018070235103", "414018070235104", "414018070235106", "414018070235108", "414018070235110", "414020070234801", "414020070240501", "414021070235201", "414024070234804", "414024070234805", "414024070234807", "414025070235201", "414032070234601")
# map all sites by type
pal <- colorFactor("Accent", domain = siteinfo$site_tp_cd)
leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(site_tp_cd),
popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = pal, title = "Site type", opacity = 1, values = ~site_tp_cd, position = "topleft")
# map all sites by treatment
pal <- colorFactor("Accent", domain = siteinfo$treatment)
leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(treatment),
popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = pal, title = "Treatment", opacity = 1, values = ~treatment, position = "topleft")
\(~\)
# parameters (nutrients samples are all filtered)
temp <- "00010" # temperature (degC)
sc <- "00095" # SpC (uS/cm)
ph <- "00400" # pH
do <- "00300" # DO (mg/L)
nh34 <- "00608" # ammonia and ammonium, as N (mg/L)
no2 <- "00613" # nitrite as N (mg/L)
no3 <- "00618" # nitrate as N (mg/L)
no23 <- "00631" # nitrate + nitrite as N (mg/L)
tn <- "62854" # total nitrogen (mg/L)
po4 <- "00671" # orthophosphate as P (mg/L)
delH <- "82082" # delta H2/H1 (per mil)
delO <- "82085" # delta O18/O16 (per mil)
pset <- c(temp, sc, ph, do, nh34, no2, no3, no23, tn, po4, delH, delO)
# RETRIEVAL
# pull water quality data
## option A (to be retired)
data <- readNWISqw(siteNumber = siteinfo$site_no, parameterCd = pset)
b <- data %>% dplyr::select(c(site_no, sample_dt, sample_tm, parm_cd, result_va))
# ## option B (suggested going forward)
# data <- readWQPqw(siteNumber = paste0("USGS-", siteinfo$site_no), parameterCd = pset)
# # select columns and rename
# b <- data %>% dplyr::select(c(MonitoringLocationIdentifier, ActivityStartDate, USGSPCode, CharacteristicName, ResultMeasure.MeasureUnitCode, ResultMeasureValue))
# names(b) <- c("site_no", "sample_dt", "parm_cd", "charname", "units", "result_va")
# b$site_no <- sub("USGS-", "", b$site_no)
# combine site info and qw data
a <- siteinfo
a$station_nm[grep("SHUBAEL", a$station_nm)] <- "SHUBAEL POND"
alldata <- right_join(a,b) %>%
mutate(param = parm_cd)
# create param names from codes
alldata$param <- dplyr::recode(alldata$param, `00010` = "temp", `00095` = "SpC", `00400` = "pH", `00300` = "DO", `00608` = "NH34", `00613` = "NO2", `00618` = "NO3", `00631` = "NO23", `62854` = "TN", `00671` = "PO4", `82082` = "delH", `82085` = "delO")
# pull water levels data
gwl <- readNWISgwl(siteNumbers = siteinfo$site_no)%>%
filter(lev_dt > "2021-09-01")
# join site info
dtw <- gwl %>%
filter(!is.na(lev_va))
dtw <- right_join(siteinfo %>% dplyr::select(station_nm, site_no, id, dec_lat_va, dec_long_va, alt_va, well_depth_va, weblink), dtw)
# SUMMARIES
# number of sample depths at each location
nwells <- alldata %>%
group_by(id) %>%
summarise(depths = length(unique(well_depth_va)))
# number of depths and samples at each location, by parameter
allsum <- alldata %>%
group_by(id, param) %>%
summarise(count = n()) %>%
pivot_wider(id_cols = id, names_from = param, values_from = count)
allsum <- left_join(allsum, nwells) %>%
dplyr::select(id, depths, temp, SpC, pH, DO, NH34:NO3, TN) %>%
arrange(depths)
# subset data for a specific sampling event
zt <- alldata %>%
dplyr::filter(sample_dt > "2021-09-01")
# estimate depth below water table (satdepth) for each sample
zt <- left_join(zt, dtw %>% select(c(site_no, lev_va))) # first, join depth to water
mls <- zt %>% filter(site_tp_cd == "multilevel sampler (MLS)") %>% # for mls, assign depth of most shallow sampled port, -0.5ft, as estimate of depth to water
group_by(id) %>%
summarise(dtw = min(well_depth_va)-0.5)
for(i in 1:nrow(mls)){zt$lev_va[zt$id==mls$id[i]] = mls$dtw[i]}
zt <- zt %>%
mutate(satdepth = well_depth_va - lev_va)
# means over shallow sample depths (<10 ft below water table) at each location, single time
zmean <- zt %>%
filter(!station_nm=="MA-A1W 517-0036") %>% # remove duplicate shallow well at A517
filter(satdepth < 10) %>%
group_by(id, parm_cd) %>%
summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 1))
# means over shallow sample depths (<10 ft below water table) at each location, time series
tsmean <- alldata %>%
filter(site_no %in% sn) %>%
group_by(id, sample_dt, param) %>%
summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 1), treatment=first(treatment))
# max over sample depths at each location
zmax <- zt %>%
group_by(id, parm_cd) %>%
summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=max(result_va, na.rm=TRUE))
# calculate percentages of N in different forms
ftn <- pivot_wider(zt, id_cols = c(id, site_no:dec_long_va, well_depth_va), names_from = param, values_from = result_va)%>%
filter(!is.na(TN)) %>%
mutate(percent_organicN = round(100*(1 - ((NH34+NO23)/TN)), 1), percent_NO23 = round(100*NO23/TN, 1))
kable(ftn %>% arrange(id) %>% dplyr::select(station_nm, site_tp_cd, temp, SpC, pH, DO, PO4, NH34, NO3, NO23, TN, percent_organicN, percent_NO23), digits = 2)
| station_nm | site_tp_cd | temp | SpC | pH | DO | PO4 | NH34 | NO3 | NO23 | TN | percent_organicN | percent_NO23 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MA-A1W 505-0036 | well cluster | 13.4 | 150 | 5.0 | 9.30 | NA | 0.02 | 4.32 | 4.32 | 4.57 | 5.0 | 94.5 |
| MA-A1W 509-0019 | water table well | 18.1 | 126 | 6.1 | 0.67 | 0 | 0.02 | 0.05 | 0.05 | 0.13 | 46.9 | 37.7 |
| MA-A1W 515-0036 | water table well | 18.8 | 77 | 5.2 | 7.80 | NA | 0.02 | 1.70 | 1.70 | 1.88 | 8.5 | 90.4 |
| MA-A1W 517-0037 | well cluster | 19.4 | 138 | 4.7 | 4.80 | NA | 0.23 | 4.02 | 4.02 | 4.58 | 7.2 | 87.8 |
| MA-A1W 521-0041 | water table well | 14.6 | 127 | 5.2 | 8.70 | NA | 0.02 | 2.16 | 2.16 | 2.40 | 9.2 | 90.0 |
| MA-A1W 522-M01-06WT | multilevel sampler (MLS) | 18.1 | 59 | 5.5 | 0.60 | 0 | 0.02 | 1.71 | 1.71 | 1.89 | 8.5 | 90.5 |
| MA-A1W 524-0039 | well cluster | 17.4 | 151 | 5.0 | 5.00 | NA | 0.02 | 3.21 | 3.21 | 3.52 | 8.2 | 91.2 |
| MA-A1W 526-0034 | well cluster | 16.5 | 116 | 5.8 | 0.69 | NA | 0.02 | 0.12 | 0.12 | 0.23 | 40.9 | 50.4 |
\(~\)
\(~\)
\(~\)
# function to make the maps
mapparam <- function(data, palette, title){
leaflet(data, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(result_va),
popup = ~paste(station_nm,"<br>","Value",result_va,"<br>","<a href='"
, weblink
, "' target='_blank'>"
, "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = palette, title = title, opacity = 1, values = ~result_va, position = "topleft")
}
z <- zmean
ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1
ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2
ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3
ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4
ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)
m5 <- mapparam(mapdata, pal, ttl)
#m5
ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6
ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis", domain = mapdata$result_va)
m7 <- mapparam(mapdata, pal, ttl)
#m7
ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)
m8 <- mapparam(mapdata, pal, ttl)
#m8
# map depth to water
pal <- colorNumeric("viridis", domain = dtw$lev_va)
m9 <- leaflet(dtw, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(lev_va),
popup = ~paste(station_nm,"<br>","Value",lev_va,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = pal, title = "Depth to water", opacity = 1, values = ~lev_va, position = "topleft")
sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)
\(~\)
\(~\)
\(~\)
z <- zmax
ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1
ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2
ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3
ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4
ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)
m5 <- mapparam(mapdata, pal, ttl)
#m5
ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6
ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis", domain = mapdata$result_va)
m7 <- mapparam(mapdata, pal, ttl)
#m7
ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)
m8 <- mapparam(mapdata, pal, ttl)
#m8
sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)
\(~\)
\(~\)
\(~\)
ttl = "Nitrate, as N, (mg/L) in wells (and pond)"
p <- ggplot(alldata %>% filter(site_tp_cd %in% c("well cluster", "water table well", "surface water (pond)")) %>% filter(parm_cd==no3), aes(x = sample_dt, y = result_va)) +
geom_point(aes(color=station_nm)) +
geom_line(aes(color=station_nm)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
ttl = "Nitrate, as N, (mg/L) in multilevel samplers"
p <- ggplot(alldata %>% filter(site_tp_cd %in% c("multilevel sampler (MLS)")) %>% filter(parm_cd==no3), aes(x = sample_dt, y = result_va)) +
geom_point(aes(color=station_nm)) +
geom_line(aes(color=station_nm)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
ttl = "Nitrate, as N, (mg/L) shallow wells, by treatment area"
p <- ggplot(alldata %>% filter(site_no %in% sn & param=="NO3"), aes(x = sample_dt, y = result_va, group = station_nm)) +
geom_point(aes(color = treatment)) +
geom_line(aes(color = treatment)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
ttl = "Nitrate, as N, (mg/L) shallow wells (mean for site), by treatment area"
p <- ggplot(tsmean %>% filter(param=="NO3"), aes(x = sample_dt, y = result_va, group = station_nm)) +
geom_point(aes(color = treatment)) +
geom_line(aes(color = treatment)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
\(~\)
\(~\)
parameters <- c( "SpC", "NO3") #"DO", "pH", "delO", "temp"
profiledata <- alldata %>%
filter(id %in% c(505, 508, 517, 522, 523, 524, 525, 526)) %>%
filter(param %in% parameters) %>%
mutate(date = as.character(sample_dt)) %>%
mutate(result_va = ifelse(param == "SpC", result_va/10, result_va)) %>%
mutate(param = replace(param, param == "SpC", "SpC/10")) %>%
arrange(id)
xr <- c(0, 50)
yr <- c(80, 0)
# function to make vertical profiles
zprofile <- function(data, title, xrange, yrange){
p <- ggplot(data %>% group_by(sample_dt), aes(x = result_va, y = well_depth_va)) +
geom_path(aes(color = param, linetype = date))+
scale_x_continuous(limits = xrange)+
theme_bw()+
scale_y_reverse(limits = yrange)+
ggtitle(title)
ggplotly(p)
}
# iterated over sites
p <- htmltools::tagList()
for (i in unique(profiledata$id)){
a <- filter(profiledata, id == i)
if(i==508){a <- filter(a, site_tp_cd == "multilevel sampler (MLS)")}
ttl <- paste0(a$site_tp_cd[1], " at ", i)
p[[i]] <- zprofile(a, ttl, xr, yr)
}
p